Sea Surface Altimetry Data Analysis¶
For this example we will use gridded sea-surface altimetry data from The Copernicus Marine Environment. This is a widely used dataset in physical oceanography and climate.
The dataset has been extracted from Copernicus and stored in google cloud storage in xarray-zarr format. It is catalogues in the Pangeo Cloud Catalog at https://catalog.pangeo.io/browse/master/ocean/sea_surface_height/
import hvplot.xarray
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
plt.rcParams['figure.figsize'] = (15, 10)
%matplotlib inline
Initialize Dataset¶
Here we load the dataset from the zarr store. Note that this very large dataset initializes nearly instantly, and we can see the full list of variables and coordinates, including metadata for each variable.
from intake import open_catalog
cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean.yaml")
ds = cat["sea_surface_height"].to_dask()
ds
<xarray.Dataset>
Dimensions: (latitude: 720, longitude: 1440, nv: 2, time: 8901)
Coordinates:
crs int32 ...
lat_bnds (time, latitude, nv) float32 dask.array<chunksize=(5, 720, 2), meta=np.ndarray>
* latitude (latitude) float32 -89.875 -89.625 -89.375 ... 89.625 89.875
lon_bnds (longitude, nv) float32 dask.array<chunksize=(1440, 2), meta=np.ndarray>
* longitude (longitude) float32 0.125 0.375 0.625 ... 359.375 359.625 359.875
* nv (nv) int32 0 1
* time (time) datetime64[ns] 1993-01-01 1993-01-02 ... 2017-05-15
Data variables:
adt (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
err (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
sla (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
ugos (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
ugosa (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
vgos (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
vgosa (time, latitude, longitude) float64 dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
Attributes:
Conventions: CF-1.6
Metadata_Conventions: Unidata Dataset Discovery v1.0
cdm_data_type: Grid
comment: Sea Surface Height measured by Altimetry...
contact: servicedesk.cmems@mercator-ocean.eu
creator_email: servicedesk.cmems@mercator-ocean.eu
creator_name: CMEMS - Sea Level Thematic Assembly Center
creator_url: http://marine.copernicus.eu
date_created: 2014-02-26T16:09:13Z
date_issued: 2014-01-06T00:00:00Z
date_modified: 2015-11-10T19:42:51Z
geospatial_lat_max: 89.875
geospatial_lat_min: -89.875
geospatial_lat_resolution: 0.25
geospatial_lat_units: degrees_north
geospatial_lon_max: 359.875
geospatial_lon_min: 0.125
geospatial_lon_resolution: 0.25
geospatial_lon_units: degrees_east
geospatial_vertical_max: 0.0
geospatial_vertical_min: 0.0
geospatial_vertical_positive: down
geospatial_vertical_resolution: point
geospatial_vertical_units: m
history: 2014-02-26T16:09:13Z: created by DUACS D...
institution: CLS, CNES
keywords: Oceans > Ocean Topography > Sea Surface ...
keywords_vocabulary: NetCDF COARDS Climate and Forecast Stand...
license: http://marine.copernicus.eu/web/27-servi...
platform: ERS-1, Topex/Poseidon
processing_level: L4
product_version: 5.0
project: COPERNICUS MARINE ENVIRONMENT MONITORING...
references: http://marine.copernicus.eu
source: Altimetry measurements
ssalto_duacs_comment: The reference mission used for the altim...
standard_name_vocabulary: NetCDF Climate and Forecast (CF) Metadat...
summary: SSALTO/DUACS Delayed-Time Level-4 sea su...
time_coverage_duration: P1D
time_coverage_end: 1993-01-01T12:00:00Z
time_coverage_resolution: P1D
time_coverage_start: 1992-12-31T12:00:00Z
title: DT merged all satellites Global Ocean Gr...xarray.Dataset
- latitude: 720
- longitude: 1440
- nv: 2
- time: 8901
- crs()int32...
- comment :
- This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system.
- grid_mapping_name :
- latitude_longitude
- inverse_flattening :
- 298.257
- semi_major_axis :
- 6378136.3
array(-2147483647, dtype=int32)
- lat_bnds(time, latitude, nv)float32dask.array<chunksize=(5, 720, 2), meta=np.ndarray>
- comment :
- latitude values at the north and south bounds of each pixel.
- units :
- degrees_north
Array Chunk Bytes 51.27 MB 28.80 kB Shape (8901, 720, 2) (5, 720, 2) Count 1782 Tasks 1781 Chunks Type float32 numpy.ndarray - latitude(latitude)float32-89.875 -89.625 ... 89.625 89.875
- axis :
- Y
- bounds :
- lat_bnds
- long_name :
- Latitude
- standard_name :
- latitude
- units :
- degrees_north
- valid_max :
- 89.875
- valid_min :
- -89.875
array([-89.875, -89.625, -89.375, ..., 89.375, 89.625, 89.875], dtype=float32) - lon_bnds(longitude, nv)float32dask.array<chunksize=(1440, 2), meta=np.ndarray>
- comment :
- longitude values at the west and east bounds of each pixel.
- units :
- degrees_east
Array Chunk Bytes 11.52 kB 11.52 kB Shape (1440, 2) (1440, 2) Count 2 Tasks 1 Chunks Type float32 numpy.ndarray - longitude(longitude)float320.125 0.375 ... 359.625 359.875
- axis :
- X
- bounds :
- lon_bnds
- long_name :
- Longitude
- standard_name :
- longitude
- units :
- degrees_east
- valid_max :
- 359.875
- valid_min :
- 0.125
array([1.25000e-01, 3.75000e-01, 6.25000e-01, ..., 3.59375e+02, 3.59625e+02, 3.59875e+02], dtype=float32) - nv(nv)int320 1
- comment :
- Vertex
- units :
- 1
array([0, 1], dtype=int32)
- time(time)datetime64[ns]1993-01-01 ... 2017-05-15
- axis :
- T
- long_name :
- Time
- standard_name :
- time
array(['1993-01-01T00:00:00.000000000', '1993-01-02T00:00:00.000000000', '1993-01-03T00:00:00.000000000', ..., '2017-05-13T00:00:00.000000000', '2017-05-14T00:00:00.000000000', '2017-05-15T00:00:00.000000000'], dtype='datetime64[ns]')
- adt(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- comment :
- The absolute dynamic topography is the sea surface height above geoid; the adt is obtained as follows: adt=sla+mdt where mdt is the mean dynamic topography; see the product user manual for details
- grid_mapping :
- crs
- long_name :
- Absolute dynamic topography
- standard_name :
- sea_surface_height_above_geoid
- units :
- m
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - err(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- comment :
- The formal mapping error represents a purely theoretical mapping error. It mainly traduces errors induced by the constellation sampling capability and consistency with the spatial/temporal scales considered, as described in Le Traon et al (1998) or Ducet et al (2000)
- grid_mapping :
- crs
- long_name :
- Formal mapping error
- units :
- m
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - sla(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- comment :
- The sea level anomaly is the sea surface height above mean sea surface; it is referenced to the [1993, 2012] period; see the product user manual for details
- grid_mapping :
- crs
- long_name :
- Sea level anomaly
- standard_name :
- sea_surface_height_above_sea_level
- units :
- m
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - ugos(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- grid_mapping :
- crs
- long_name :
- Absolute geostrophic velocity: zonal component
- standard_name :
- surface_geostrophic_eastward_sea_water_velocity
- units :
- m/s
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - ugosa(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- comment :
- The geostrophic velocity anomalies are referenced to the [1993, 2012] period
- grid_mapping :
- crs
- long_name :
- Geostrophic velocity anomalies: zonal component
- standard_name :
- surface_geostrophic_eastward_sea_water_velocity_assuming_sea_level_for_geoid
- units :
- m/s
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - vgos(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- grid_mapping :
- crs
- long_name :
- Absolute geostrophic velocity: meridian component
- standard_name :
- surface_geostrophic_northward_sea_water_velocity
- units :
- m/s
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray - vgosa(time, latitude, longitude)float64dask.array<chunksize=(5, 720, 1440), meta=np.ndarray>
- comment :
- The geostrophic velocity anomalies are referenced to the [1993, 2012] period
- grid_mapping :
- crs
- long_name :
- Geostrophic velocity anomalies: meridian component
- standard_name :
- surface_geostrophic_northward_sea_water_velocity_assuming_sea_level_for_geoid
- units :
- m/s
Array Chunk Bytes 73.83 GB 41.47 MB Shape (8901, 720, 1440) (5, 720, 1440) Count 1782 Tasks 1781 Chunks Type float64 numpy.ndarray
- Conventions :
- CF-1.6
- Metadata_Conventions :
- Unidata Dataset Discovery v1.0
- cdm_data_type :
- Grid
- comment :
- Sea Surface Height measured by Altimetry and derived variables
- contact :
- servicedesk.cmems@mercator-ocean.eu
- creator_email :
- servicedesk.cmems@mercator-ocean.eu
- creator_name :
- CMEMS - Sea Level Thematic Assembly Center
- creator_url :
- http://marine.copernicus.eu
- date_created :
- 2014-02-26T16:09:13Z
- date_issued :
- 2014-01-06T00:00:00Z
- date_modified :
- 2015-11-10T19:42:51Z
- geospatial_lat_max :
- 89.875
- geospatial_lat_min :
- -89.875
- geospatial_lat_resolution :
- 0.25
- geospatial_lat_units :
- degrees_north
- geospatial_lon_max :
- 359.875
- geospatial_lon_min :
- 0.125
- geospatial_lon_resolution :
- 0.25
- geospatial_lon_units :
- degrees_east
- geospatial_vertical_max :
- 0.0
- geospatial_vertical_min :
- 0.0
- geospatial_vertical_positive :
- down
- geospatial_vertical_resolution :
- point
- geospatial_vertical_units :
- m
- history :
- 2014-02-26T16:09:13Z: created by DUACS DT - 2015-11-10T19:42:51Z: Change of some attributes - 2017-01-06 12:12:12Z: New format for CMEMSv3
- institution :
- CLS, CNES
- keywords :
- Oceans > Ocean Topography > Sea Surface Height
- keywords_vocabulary :
- NetCDF COARDS Climate and Forecast Standard Names
- license :
- http://marine.copernicus.eu/web/27-service-commitments-and-licence.php
- platform :
- ERS-1, Topex/Poseidon
- processing_level :
- L4
- product_version :
- 5.0
- project :
- COPERNICUS MARINE ENVIRONMENT MONITORING SERVICE (CMEMS)
- references :
- http://marine.copernicus.eu
- source :
- Altimetry measurements
- ssalto_duacs_comment :
- The reference mission used for the altimeter inter-calibration processing is Topex/Poseidon between 1993-01-01 and 2002-04-23, Jason-1 between 2002-04-24 and 2008-10-18, OSTM/Jason-2 since 2008-10-19.
- standard_name_vocabulary :
- NetCDF Climate and Forecast (CF) Metadata Convention Standard Name Table v37
- summary :
- SSALTO/DUACS Delayed-Time Level-4 sea surface height and derived variables measured by multi-satellite altimetry observations over Global Ocean.
- time_coverage_duration :
- P1D
- time_coverage_end :
- 1993-01-01T12:00:00Z
- time_coverage_resolution :
- P1D
- time_coverage_start :
- 1992-12-31T12:00:00Z
- title :
- DT merged all satellites Global Ocean Gridded SSALTO/DUACS Sea Surface Height L4 product and derived variables
Visually Examine Some of the Data¶
Let’s do a sanity check that the data looks reasonable. Here we use the hvplot interactive plotting library.
ds.sla.hvplot.image(
'longitude',
'latitude',
rasterize=True,
dynamic=True,
width=800,
height=450,
widget_type='scrubber',
widget_location='bottom',
cmap='RdBu_r',
)
Create and Connect to Dask Distributed Cluster¶
from dask.distributed import Client
from dask_gateway import Gateway
gateway = Gateway()
cluster = gateway.new_cluster()
cluster.adapt(minimum=1, maximum=20)
cluster
** ☝️ Don’t forget to click the link above to view the scheduler dashboard! **
client = Client(cluster)
client
Client
|
Cluster
|
Timeseries of Global Mean Sea Level¶
Here we make a simple yet fundamental calculation: the rate of increase of global mean sea level over the observational period.
# the number of GB involved in the reduction
ds.sla.nbytes / 1e9
73.8284544
# the computationally intensive step
sla_timeseries = ds.sla.mean(dim=('latitude', 'longitude')).load()
sla_timeseries.plot(label='full data')
sla_timeseries.rolling(time=365, center=True).mean().plot(label='rolling annual mean')
plt.ylabel('Sea Level Anomaly [m]')
plt.title('Global Mean Sea Level')
plt.legend()
plt.grid()
In order to understand how the sea level rise is distributed in latitude, we can make a sort of Hovmöller diagram.
sla_hov = ds.sla.mean(dim='longitude').load()
fig, ax = plt.subplots(figsize=(12, 4))
sla_hov.name = 'Sea Level Anomaly [m]'
sla_hov.transpose().plot(vmax=0.2, ax=ax)
<matplotlib.collections.QuadMesh at 0x7f8052a09ed0>
We can see that most sea level rise is actually in the Southern Hemisphere.
Sea Level Variability¶
We can examine the natural variability in sea level by looking at its standard deviation in time.
sla_std = ds.sla.std(dim='time').load()
sla_std.name = 'Sea Level Variability [m]'
ax = sla_std.plot()
_ = plt.title('Sea Level Variability')
